1. Basic concept of programming, data types, vectors and factors

2. Matrices and data.frames, data import/export, graphics and visualization

3. Tidy data, data wrangling

4. Today: Spatial data in R

library(raster)

ras <- raster(ncol = 100, nrow = 100, xmn = 0, xmx = 100, ymn = 0, ymx = 100)

ras

## class       : RasterLayer 
## dimensions  : 100, 100, 10000  (nrow, ncol, ncell)
## resolution  : 1, 1  (x, y)
## extent      : 0, 100, 0, 100  (xmin, xmax, ymin, ymax)
## coord. ref. : NA

m <- rnorm(n = ncell(ras), mean = 0, sd = 1)

values(ras) <- m

plot(ras)

ras2 <- ras + 10

ras3 <- ras * ras2

ras4 <- ras^2

ras5 <- (ras - cellStats(ras, stat = mean)) / 
                                cellStats(ras, stat = sd)

ras.na <- ras

ras.na[ras.na > 0.5] <- NA

k <- matrix(c(  0, 0.5, 1,
              0.5, 1.5, 2,
              1.5, 6.0, 3), 
            ncol = 3, byrow = TRUE)
k

##      [,1] [,2] [,3]
## [1,]  0.0  0.5    1
## [2,]  0.5  1.5    2
## [3,]  1.5  6.0    3

ras.reclass <- reclassify(x = ras, rcl = k)

ras[10, 20]

##           
## 0.7534085

ras[10, 20] <- 10

ras.stack <- stack(ras, ras2)

ras.stack

## class       : RasterStack 
## dimensions  : 100, 100, 10000, 2  (nrow, ncol, ncell, nlayers)
## resolution  : 1, 1  (x, y)
## extent      : 0, 100, 0, 100  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## names       :   layer.1,   layer.2 
## min values  : -3.592167,  6.407833 
## max values  :  10.00000,  13.60653

layer1 <- subset(x = ras.stack, subset = 1)

layer1

## class       : RasterLayer 
## dimensions  : 100, 100, 10000  (nrow, ncol, ncell)
## resolution  : 1, 1  (x, y)
## extent      : 0, 100, 0, 100  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## data source : in memory
## names       : layer.1 
## values      : -3.592167, 10  (min, max)

data <- as.data.frame(ras.stack)

head(data)

##        layer.1   layer.2
## 1 -0.531517486  9.468483
## 2 -0.183922864  9.816077
## 3  0.003593095 10.003593
## 4 -0.989860242  9.010140
## 5 -3.592166533  6.407833
## 6 -1.795983887  8.204016

elevation <- raster("Data/elevation.envi")

elevation

## class       : RasterLayer 
## dimensions  : 709, 1306, 925954  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : 16.90833, 27.79167, 44.33333, 50.24167  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /Users/corneliussenf/Dropbox/Teching/Quantitative Methods/Week-13/Data/elevation.envi 
## names       : elevation

projection(elevation)

## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

plot(elevation)

# use stack() or brick() for multi-layer raster files
climate <- stack("Data/climate.envi")

# get the number of raster layers
nlayers(climate)

## [1] 4

# get the raster layer names
names(climate)

## [1] "X.Tmax."  "X.Tmin."  "X.Tmean." "X.Prec."

# rename the raster layers
names(climate) <- c("Tmax", "Tmin", "Tmean", "Prec")

writeRaster(elevation, 
            filename= "Data/elev.envi",
            overwrite = TRUE)

library(sp)

coords_srs <- data.frame(x = runif(100, 17, 27), 
                         y = runif(100, 45, 50))

pts_srs <- SpatialPoints(coords = coords_srs)
pts_srs

## class       : SpatialPoints 
## features    : 100 
## extent      : 17.04253, 26.99759, 45.01899, 49.98391  (xmin, xmax, ymin, ymax)
## coord. ref. : NA

proj <- projection(elevation)

CRS(proj)

## CRS arguments:
##  +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

pts_srs <- SpatialPoints(coords = coords_srs,
                    proj4string = CRS(proj))

plot(elevation)
plot(pts_srs, add = TRUE)

library(mapview)

mapview(climate)

df.climate <- extract(climate, pts_srs, df = TRUE)

class(df.climate)

## [1] "data.frame"

head(df.climate)

##   ID     Tmax     Tmin     Tmean     Prec
## 1  1 161.3333 60.66667 110.66666 47.83333
## 2  2 137.1667 38.58333  87.58334 47.33333
## 3  3 118.1667 38.00000  77.83334 56.33333
## 4  4 154.0000 53.91667 103.58334 46.75000
## 5  5 116.8333 28.33333  72.41666 65.33334
## 6  6 147.9167 40.50000  93.91666 48.50000

spdf.climate <- extract(climate, pts_srs, sp = TRUE)

spdf.climate

## class       : SpatialPointsDataFrame 
## features    : 100 
## extent      : 17.04253, 26.99759, 45.01899, 49.98391  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 4
## names       :            Tmax,             Tmin,            Tmean,             Prec 
## min values  :           83.75,             -2.5, 40.3333320617676, 43.4166679382324 
## max values  : 164.16667175293, 65.6666641235352,           114.25, 86.4166641235352

You can acess SpatialPointsDataFrame the same way as data.frame!

head(spdf.climate)

##       Tmax     Tmin     Tmean     Prec
## 1 161.3333 60.66667 110.66666 47.83333
## 2 137.1667 38.58333  87.58334 47.33333
## 3 118.1667 38.00000  77.83334 56.33333
## 4 154.0000 53.91667 103.58334 46.75000
## 5 116.8333 28.33333  72.41666 65.33334
## 6 147.9167 40.50000  93.91666 48.50000

spdf.climate$Tmax[1]

## [1] 161.3333

library(rgdal)

pts.sample <- readOGR('Data', 'pts_sample')

## OGR data source with driver: ESRI Shapefile 
## Source: "Data", layer: "pts_sample"
## with 142 features
## It has 1 fields
## Integer64 fields read as strings:  id

pts.sample

## class       : SpatialPointsDataFrame 
## features    : 142 
## extent      : 17.80209, 27.04407, 44.64419, 49.84834  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 1
## names       : id 
## min values  :  0 
## max values  : 99

writeOGR(spdf.climate, "Data", "pts_climate", driver = "ESRI Shapefile", overwrite_layer = TRUE)